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Abstract 


The Refined Zigzag Theory (RZT) for homogeneous, laminated 
composite, and sandwich plates is revisited to offer a fresh insight into 
its fundamental assumptions and practical possibilities. The theory is 
introduced from a multi-scale formalism starting with the inplane 
displacement field expressed as a superposition of coarse and fine 
contributions. The coarse displacement field is that of first-order shear- 
deformation theory, whereas the fine displacement field has a piecewise- 
linear zigzag distribution through the thickness. The resulting kinematic 
field provides a more realistic representation of the deformation states of 
transverse-shear-flexible plates than other similar theories. The 
condition of limiting homogeneity of transverse-shear properties is 
proposed and yields four distinct sets of zigzag functions. Analytic 
solutions for highly heterogeneous sandwich plates undergoing 
elastostatic deformations are used to identify the best-performing zigzag 
functions. Unlike previously used methods, which often result in 
anomalous conditions and nonphysical solutions, the present theory does 
not rely on transverse-shear-stress equilibrium constraints. For all 
material systems, there are no requirements for use of transverse-shear 
correction factors to yield accurate results. To model homogeneous 
plates with the full power of zigzag kinematics, infinitesimally small 
perturbations in the transverse shear properties are derived, thus 
enabling highly accurate predictions of homogeneous-plate behavior 
without the use of shear correction factors. The RZT predictive 
capabilities to model highly heterogeneous sandwich plates are 
critically assessed, demonstrating its superior efficiency, accuracy, and 
a wide range of applicability. The present theory, which is derived from 
the virtual work principle, is well-suited for developing computationally 
efficient, C° -continuous finite elements, and is thus appropriate for the 
analysis and design of high-performance load-bearing aerospace 
structures. 


Nomenclature 

a 

c° 

c u and Co- 


lateral dimension of a square plate 

a continuous function of (x,,x 2 ) coordinates whose first-order 
derivatives are discontinuous along finite element interfaces 
intersections of the cylindrical edge surfaces ( S u , S a ) with the 
middle surface S m where displacements and traction resultants 
are prescribed, respectively 
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a continuous function of z coordinate whose first-order 
derivative is discontinuous along material-layer interfaces 
inplane and transverse shear elastic stiffness coefficients for the 
A*th layer 

Y oung's moduli of the Arth layer 

transverse shear piecewise constant functions 

weighted-average, laminate-dependent transverse shear constants 

shear moduli of the Alh layer 

total plate (laminate) thickness 

thickness of the /cth material layer 

shear correction factor for First-order Shear Deformation Theory 
(FSDT) 

number of material layers through the plate thickness 
applied transverse pressure [force/length 2 ] 

reference middle plane of the laminate 

parts of the cylindrical edge surface of the laminate where 

displacements and tractions are prescribed, respectively 

transverse-shear compliance coefficients 

unit outward tangential and normal vectors to the mid-plane 

boundary 

prescribed inplane and transverse shear tractions 

inplane and transverse components of the displacement vector in 

the /cth material layer 

coarse kinematic description 

fine kinematic description 

central transverse displacement averaged across plate thickness 
normalized inplane displacement, normal stress, and transverse 
shear stress 

kinematic variables of the refined zigzag plate theory 

reference coordinates with axes positioned in the middle plane of 

plate 

thickness coordinate axis 
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thickness coordinate of the interface between the Mi and (& + 1) 
layers 

subscript (« = 1,2) 
average transverse shear strains 

transverse shear strain function corresponding to cylindrical 
bending problem 
variational operator 
partial differentiation 

small, dimensionless piecewise constant function 
strains in the Mi layer 

dimensionless thickness coordinates of the Mi material layer 
piecewise constant bending rotation 

the angle between the principal material direction x[ and the 

x\ axis, formed by the right-hand rule 

dimensionless piecewise constant function 

Poisson ratios of the Mi layer 

stresses in the Alh layer 

zigzag functions in the Mi layer 

dimensionless inplane displacements along the interface between 

the Mi and (k + 1) layers 

amplitude functions of zigzag displacements 
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1. Introduction 


Advances in composites technology over the past four decades spurred an ever increasing use of 
composite materials in civilian and military aircraft, aerospace vehicles, naval and civil structures. 
Application of composite materials leads to lightweight structures that offer high-performance and long- 
term durability and reliability. The newest state-of-the-art civilian aircraft, the Boeing 787, is fifty percent 
composite and uses carbon-epoxy materials for its fuselage, wings, and many load-bearing components. 

Stress-analysis methods based on classical assumptions, which suppress the transverse shear and 
normal effects, have significant limitations particularly when applied to multilayered composites 
undergoing bending and transverse shear deformations. Load-bearing thick-section composites and 
sandwich structures can exhibit significant deformations due to transverse shear and thickness-stretch 
effects, especially in regions of stress concentration and when undergoing high-frequency responses. 
Many structural theories of first and higher order and their finite element implementations have been 
developed. The most widely used finite elements for plate and shell bending are those based on First- 
order Shear-Deformation Theory (FSDT) [1-4], because of their computational efficiency and relatively 
wide range of applicability. Nevertheless, it is well recognized that such models have the tendency to 
underestimate the inplane strains and stresses, particularly in highly heterogeneous and thick composite 
and sandwich laminates, and overestimate natural frequencies of relatively high-frequency vibration 
modes [5-8]. Moreover, FSDT’s accuracy is related to the use of shear correction factors that are material 
lay-up dependent. Higher-order theories, that take into account transverse shear deformability [8], are 
more accurate than FSDT and normally do not require shear correction factors. Theories based on the 
equivalent-single-layer assumption [9-12] also offer reduced computational complexity; however, they 
fail to model the zigzag-shaped cross-sectional distortion typical of heterogeneous laminates. Layer-wise 
theories [13, 14] have quasi-three-dimensional predictive capabilities; however, the computational effort 
is excessively great for most practical applications. To model sandwich laminates with sufficient accuracy 
commonly requires application of special sandwich theories to account for a zigzag-shaped, through-the- 
thickness distribution of inplane displacements and strains, at the expense of a much higher modeling and 
kinematic complexity [15-17]. Alternatively, three-dimensional finite element models are used to analyze 
relatively small structural regions of interest, and these are often coupled with two-dimensional 
discretizations to achieve smaller and computationally feasible models. Computational methods for 
nonlinear progressive failure analysis commonly employ cohesive elements and ply-by-ply 
discretizations, resulting in exceedingly intensive computations, particularly when large structures are 
analyzed [18, 19]. 

Zigzag theories pioneered by Di Sciuva, both in linear [20-23] and in cubic versions [24, 25], may be 
considered as a good compromise between adequate accuracy and low computational cost. The key idea 

is to add a piecewise-linear, zigzag-shaped (i.e., C? -continuous) contribution to a globally linear or cubic 
through-the-thickness distribution of the inplane displacements. The superposed zigzag kinematics are 
determined in such a way as to satisfy equilibrium of the transverse shear stresses through the laminate 
thickness, as required by the theory of elasticity. The formulation results in a fixed number of kinematic 
unknowns, equal to that of FSDT, which does not depend on the number of layers. These theories often 
yield response predictions comparable to those of layer-wise theories; however, various flaws inherent in 
these theories have prevented their acceptance in practical applications. 
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Averill [26] recognized two major drawbacks that plague many previously mentioned zigzag theories: 
(i) C 1 -continuous functions are required to approximate the deflection variable within the finite element 
framework - the type of approximations that are especially undesirable for plate and shell finite elements, 
and (ii) transverse shear stresses calculated from constitutive equations vanish erroneously along clamped 
boundaries. 

Building on the observations of Averill [26], Tessler et. al. presented [27-31] a Refined Zigzag Theory 
(RZT) that augments FSDT with a novel zigzag representation of the inplane displacements. The 
kinematic field does not require enforcement of transverse-shear-stress continuity to yield accurate 
results. Both drawbacks of the original zigzag theories are overcome because: (i) only first derivatives of 

the kinematic unknowns are present in the definition of the strain field, thus leading to C° -continuous 
shape functions for beam and plate finite elements [32-35], and (ii) all field equations (i.e., equilibrium 
equations, constitutive equations, boundary conditions, and strain-displacement relations) are consistently 
derived from the virtual work principle without engendering any transverse-shear-force anomalies. Since 
the transverse shear forces are fully consistent with respect to the physical and variational requirements, 
they do not vanish erroneously along clamped boundaries. Furthermore, the theory has been demonstrated 
to yield consistently superior results over a wide range of aspect ratios and material systems, including 
thick laminates with a high degree of transverse shear flexibility, anisotropy, and heterogeneity. 

In this paper, the Refined Zigzag Theory [30] is revisited to reexamine its fundamental assumptions 
and to explore new modeling possibilities. In Section 2, the theory is introduced from a multi-scale 
formalism. The kinematic field is represented as a superposition of coarse and fine kinematic 
descriptions , where the coarse kinematics are those of FSDT and the fine kinematics are represented by 

piecewise C? -continuous zigzag functions through the thickness. The definition of the zigzag functions is 
further refined by assuming that the coarse kinematics form the average laminate response only, whereas 
the fine kinematics represent a perturbation from the average response. The strain-displacement equations 
and constitutive relations are presented in Section 3. In Section 4, the zigzag functions are defined in 
terms of layer-level thickness coordinates. The methodology for determining the zigzag functions from 
the transverse-shear constitutive relations is reviewed following the approach in [30]. Then, a new 
perspective, based on the concept of limiting homogeneity , is proposed for the determination of zigzag 
functions. This new methodology permits four distinct sets of zigzag functions to be obtained, one of 
which is coincident with the form derived in [30]. The concept of limiting homogeneity is further invoked 
for the purpose of modeling homogeneous plates with the full power of zigzag kinematics. Analytic 
results for homogeneous orthotropic plates and highly heterogeneous sandwich laminates are presented in 
Section 5 to demonstrate RZT’s wide range of applicability. 


2. Kinematics description 


Consider a multilayered composite plate of uniform thickness 2 h composed of perfectly bonded 
orthotropic layers (or laminae) as shown in Figure 1. Points of the plate are located by the orthogonal 
Cartesian coordinates (Ay,x ? ,z). The ordered pair (ij,x 2 ) e S m denotes the inplane coordinates, where 

S m represents the set of points given by the intersection of the plate with the plane z = 0 , referred to 
herein as the middle reference plane (or midplane). The symbol z e [-/z, h] denotes the domain of the 
through-the-thickness coordinate, with z = 0 identifying the plate’s midplane. The plate is subjected to a 
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normal-pressure loading, q(x\,x 2 ) , applied at the midplane, S m , that is defined as positive when acting 
in the positive z direction. In addition, a traction vector, {T x ,T 2i T , ) , is prescribed on cz S , where S 

denotes the total cylindrical-edge surface. On the remaining pail of the edge surface, S n cz S , 
displacement restraints are imposed (or prescribed). The sections of the plate edge are related by 
■S' <jS„ = S and S„ n S, , = 0 . Moreover, the curves C„ = S„ r\ S m and C„ = S„ rv S,„ define the two 

G U G ll 5 g g m u u m 

parts of the total perimeter C = C a ^jC u surrounding the midplane region, S m . Finally, it is presumed 
that small-strain assumptions are valid, and that body and inertial forces are negligible. 

The deformation of a plate subjected to applied loading and boundary restraints can be described by 
the displacement vector in terms of its Cartesian components as u = . Using a multi-scale 

formalism, the material point at (x 1? x 2 ,z) situated within the Ath material layer — which could be a ply 
or a fraction of a ply thickness, if the ply is further discretized through the thickness — is undergoing the 
displacement that may be expressed as a supeiposition of a coarse kinematic description and a fine 

kinematic description (/ = l,2,z) , i.e., 

} ( X \ ’ X 2’ z ) = ll a(c ) (*1 > *2 ’ Z ) + U a(f) (*1 > *2 , Z) (a = 1, 2) 

(1) 

u[ k) (x l9 x 2 ,z) = u z(c) (x l ,x 2 ,z) + u^ f) (*! ,x 2 ,z) 

where the superscript (k) (k = 1,..., N , where N denotes the total number of material layers through the 
plate thickness) indicates that a quantity is dependent upon the Alh layer constitutive properties; whereas, 
in the ensuing discussion, the subscript ( k ) defines quantities corresponding to the interface between the 

k and (£+1) layers. The Ath layer thickness is defined in the range [z^_^ , z (/t) ] (see Figure 1). 

To describe the coarse kinematics, first-order shear deformation theory (FSDT) can be used, in which 
case 

«1 (C) (*1 > * 2 , z) = M (*1 t x 2 ) + Z e x (X x , X 2 ) 

U 2{c) ( *1 ,X 2 ,Z) = V(X 1 ,X 2 ) + Z 0 2 (X, , X 2 ) 

w_.(x 1 ,x 2 ,z) = w(x l5 x 2 ) 

where a and v represent the uniform through-the-thickness displacements in the Xj and x 2 coordinate 
directions, respectively; 0 X and 0 2 represent average bending rotations of the transverse normal about the 
positive x 2 and the negative Xj directions, respectively, and they contribute to the linear distributions 
through the thickness; and w is the coarse or average description of transverse deflection. 

The fine kinematic description with the resolution on the scale of the Ath material layer is represented 

(k) 

by the piecewise-continuous functions u a (f) ■ general, the addition of the fine kinematic description 
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adds 3x(jV + l) independent variables. If this generality of the fine description is maintained, then Eq. 
(1) would represent what is commonly known as a layer- wise theory, where the number of unknown 
kinematic variables (or modes) is proportional to the number of material or mathematical layers that are 
modeled within such a theory. One way to simplify the fine kinematic description is first to represent the 
fine kinematics in terms of separable functions of the inplane (jq , x 2 ) and transverse (z) coordinates, 

and then impose a set of interfacial constraint conditions, with this general strategy leading to what is 
commonly known as a zigzag formulation (e.g., see Di Sciuva [20]). For many years it has been widely 
accepted that such interface constraint conditions should necessarily be conditions of perfect equilibrium 
of the transverse shear stresses along the layer interfaces. Whereas from the mechanics perspective the 
stress equilibrium needs to be perfectly satisfied, the limited kinematic freedom associated with the 
relatively simple displacements assumed a priori often lead to an over constraining of the kinematic field 
and subsequent pathological conditions and anomalies within the resulting theory. 

In Tessler et al. [30], the fine displacement description is assumed as 

“«(/) (*1 ■ * 2 ’ Z ) S O ) Va (*1 . *2 ) 

"!(/)(* i’X 2 ,z) = 0 (3) 

allowing only for the inplane displacements to have a fine description, where $„ k \z) are piecewise- 
continuous through-the-thickness zigzag functions that stem from laminate heterogeneity, and y/ a (Xj , x 2 ) 
are independent functions that may be interpreted as the amplitudes of the zigzag displacements. In its 
simplest form, as proposed by Tessler et al. [30], (f^\z) are piecewise linear, C° -continuous functions 

through the thickness, having discontinuous thickness-derivatives (z) at the layer interfaces. 

The refined zigzag theory of Tessler et al. [30] has seven kinematic variables (modes) regardless of the 
number of material layers through the thickness, with the kinematic field described by the orthogonal 
components of the displacement vector as 


u \ k \x\,x 2 , z ) = u(x l ,x 2 ) + z 0fx l ,x 2 ) + <fi k \z) i// l (x l ,x 2 ) 
u i 2 k \x l ,x 2 ,z) = v(x l ,x 2 ) + z <9 2 (Xi,x 2 ) + $ k) (z) i// 2 {x lt x 2 ) 
m 2 (x 1 ,x 2 ,z) = w(x 1 ,x 2 ) (4) 

During applied loading the material fiber, which prior to deformation is straight and perpendicular to the 
reference midplane, deforms and rotates about the x a axes. The rotations, henceforth denoted as , 
are given as 
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(a) (b) 

Figure 1. (a) General plate notation, and (b) lamination notation. 


U «l =0 a{ X l’ X 2) + ^z{ Z ) ( 5 ) 

where, henceforward, (•) a = 3(*) / dx a denotes a partial derivative with respect to the midplane 

coordinate, x a . In Eqs. (5), are piecewise constant, and @ a (x l ,x 2 ) an >d V a (*i>* 2 ) are uniform 
with respect to the z coordinate; therefore, each material layer rotates in a piecewise constant manner 
described by the rotations . Integrating Eqs. (5) over the laminate thickness and dividing by the total 
thickness gives the average rotations 



If the average rotations about the x\ and x 2 axes of the initially normal fiber are represented by the 
coarse rotations 6 a (^,* 2 ) on ^ ( as * n Di Sciuva et al. [31]), then the following integral must vanish 
identically 

l_VS( z ) rfz=o 

Integration of Eq. (7) results in the requirement that the top (surface) and bottom values of the respective 
zigzag functions are equal, i.e., 

j h _ h <tal { Z ) dz = <t { a N) { h )-<f > a ) (- h ) = 0 ( 7 - *) 
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A special case of Eq. (7.1) and a particularly convenient choice is to select the top and bottom values of 
the zigzag functions to vanish identically 


^\-h) = €\h) = 0 (8) 

The above homogeneity conditions of the zigzag functions also imply that surface displacements are only 
defined in terms of the coarse variables, i.e., 

u[ k) (x l ,x 2 ,±h) = u (Xj , x 2 ) ± h e x (x l , x 2 ) 

u { 2 k) (x l ,x 2 ,±h) = v(Xj , x 2 ) ± h e 2 (Xj , x 2 ) r9 . 


3. Strain-displacement and constitutive relations 

The inplane and transverse shear strains derived from linear strain-displacement relations are 


^n )=w j + z ^U+^ (A Vi,i (10.1) 

S 22 = V 2 + z @2,2 + 4 ^2.2 (10.2) 

Yu = U ,2 + V .l + Z (&L2 + # 2 .l) + 4 '' Vl.2 + 4 k V 2 .l ( 10 - 3 ) 

Yaz =™a + 0a + <f f alv' a (« = U) (10.4) 

The zigzag-function property given by Eq. (7) ensures that the average transverse shear strains are those 
which correspond to FSDT (a coarse description), i.e., 

r a =j - f_ h TaJdz =w a +e a (ii) 


The generalized Hooke's law for the Ath orthotropic material layer, whose principal material directions are 
arbitrary with respect to the reference coordinates (x^x-,) e S m , is expressed as 


^11 

(k) 

A 

02 

0 6 

0 

0 ■ 

(■ k ) 

^11 

ct 22 


Q 2 

(-22 

06 

0 

0 


S 22 

T n 

= 

q 6 

^26 

0* 

0 

0 

* 

Y 12 

T 2z 


0 

0 

0 

Q22 

Qn 


Y 2z 



0 

0 

0 

Qn 

0: J 


Yiz. 


( 12 ) 


9 



where C\ k } and Q ■ J k) are the transformed elastic-stiffness coefficients referred to the (x l ,x 2 ,z) 
coordinate system, relative to the stress condition that ignores the transverse normal stress. To aid 
in subsequent discussions, expressions for Q-J' ’ are given in terms of their principal transverse-shear 

moduli, (jj 3 ^ and , and an angle (the angle between the principal material direction x[ and the 
Xj axis, formed by the right-hand rule) 


Q[\ ] = cos 2 (0 m )G{P + sin 2 (0 {k) ) G[ 3 } ; Q[ k) = sin(# (Ar) ) cos(^ ) ) (Gjf ) -Gg } ); 
Q 22 = cos 2 (0 ( ‘ k ' i )G^\ ) + sin 2 (0M)G«\ 


( 12 . 1 ) 


4. Zigzag functions and transverse shear constitutive relations 

The refined zigzag functions (or zigzag displacements) of the present theory are defined by piecewise 
linear, C? -continuous functions (i.e., C° -continuous functions through the laminate thickness.) For 

convenience, the zigzag functions <f> { a k \z) ( a = 1,2) , which have units of length, are defined in terms of 
their respective layer-interface values <j> a ^ (/ = 0,1,..., N) (see Figure 2 depicting the notation for a 
three-layered laminate). 



Figure 2. Notation for a three-layered laminate and ^ k \z) zigzag functions defined in terms 
of interface displacements (f> a (^ ■ 
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(13) 


For the Ath material layer located in the range | J , the zigzag functions are given as 

4' ( «=[(z-z (i _ 1) )/A w -l]e[-U] (i = l,...,JV) (14) 

with the first layer beginning at z^ = —h , the last N th layer ending at z^=h , and the k th layer 
ending at z^ = z^_^ + 2h {k \ where 2h {k 1 denotes the k th layer thickness. 

Evaluating Eqs. (13) at the layer interfaces gives rise to the definitions of the interface displacements 

u- n =e{? k) =-i), a?) 

where, according to Eq. (8) , the top and bottom interfacial displacements are set to vanish 

< l , a(0) =( t , a(N)-Q ( 16 ) 

Differentiating Eqs. (13) with respect to the z coordinate yields the piecewise-constant functions 

& 1 = 2fo(k) ( ^»<*> “ ) ( 1 7 ) 

Straightforward algebraic manipulations provide the alternative forms for ^ A) (z) in terms of (z) 

l)) + «W-0 OS) 


4.1 Zigzag function selection according to Tessler et al. [30] 

The approach in [30] begins by expressing the transverse shear stresses in the following form 


T \z 

T 2z. 
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Qn Qn 

012 022 
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(19.1) 


or alternatively as 
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(19.2) 



In Eqs. (19.2), the second and third stress vectors include, as their normalization (or scaling) factors, 
coefficients |l + (p^l j ( a =1,2) involving their respective zigzag functions. If the normalization 
factors are selected to be constant for all material layers, i.e., 

fiS(l + C) sG a (2°) 


where G a are yet undetermined constants, then solving Eqs. (20) for (f>^\ results in 



G /O ik) - 1 

a z-aa 


( 21 ) 


Substituting Eqs. (21) into Eqs. (7) reveals that the G a constants represent some weighted-average 
transverse-shear stiffness coefficients of a laminate, given by 


( 1 r'' dz ^ 

"'_r i * 

W~ h Q%> 

UwesJ 


(22) 


Using Eqs. (17), all interfacial displacements are readily determined in terms of (ft^l as 

&(*>=&<*-., + 2 (23.1) 

Accounting for the homogeneous boundary conditions, Eqs. (16), the interfacial displacements are further 
simplified as 

«W)=Z 2/,< °C (23-2) 

i=l 

With (f)^] and G a fully defined by Eqs. (21) and (22), the definitions of the zigzag functions given by 

Eqs. (18) are complete. Invoking Eqs. (10.4) and (22), piecewise-constant distributions of the transverse 
shear strains are maintained and are given as 
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(24) 


y (k) = vv + 6 + 

/ az ,a a 


'eSf A w 

H0 


^ ^ /O(^) 
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ace y 


Note that the present methodology for deriving the zigzag functions is purely constitutive-based; in 
contrast to the previous zigzag models, there exist no kinematic constraints. In [30], extensive analytical 
studies for simply supported and clamped laminated composite and sandwich plates have been carried out 
and have demonstrated superior modeling capabilities of this theory. 

4.2 Homogeneous limit methodology for zigzag function selection 

In what follows a new methodology for obtaining suitable zigzag functions is examined. The approach 
exploits an intrinsic property of zigzag functions referred to herein as the homogeneous limit property. It 

is postulated that a zigzag function (f^ [a = 1,2; k = 1,..., Af) has a non-vanishing distribution 

provided that the layer-wise transverse shear properties on the x a cross-sections are different, i.e., when 

the material layers are heterogeneous with respect to their transverse shear properties. The converse of 
this property is that under the homogeneous limit defined as 

Qaa -» G a (25) 


where G a are the constant transverse-shear moduli corresponding to the x a planes (i.e., the planes 

normal to the x a coordinate axes), the zigzag functions are required to vanish in a limiting sense under 
all admissible deformations, i.e., 

*?> (*)->» if Ql k l^G a (26) 

The homogeneous limit conditions, Eq. (25), also imply that all material layers have the same angular 
orientation, 6^ k \ It immediately follows from Eq. (26) that </>#] also vanish in a limiting sense (refer to 
Eqs. (18)), i.e., 

f Qaa G a (27) 

Moreover, the limiting conditions described by Eqs. (26) and (27) guarantee that all interior interfacial 
displacements of the zigzag functions must vanish in a limiting sense as the laminate approaches its 
homogeneous limit 

ta ( /)->° 0 = (28) 

whereas the exterior (bottom and top) interfacial displacements have been set to zero explicitly in Eq. 
(16). Rewriting Eq. (25) by moving G a to the left-hand side of the equation, and then dividing by 

or by G a , gives two sets of dimensionless functions labeled as g ( ^ 1 and gj k J : 
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(29) 


{gi k) , g's^Gjosa- 1 , 

where, for expediency, G a has been renamed as G Ia in the definition of gfj . It can be seen that g {k) 
and gj k J have the same characteristics as (f) {k \ ; specifically, (1) these functions are dimensionless, 
piecewise-constant and dependent upon the Ath layer’s transverse shear stiffnesses, and (2) each of these 
functions vanishes in a limiting sense under the homogeneous limit, Eq. (25). Therefore, either g (k) and 

gj k J can be considered to represent (f (k \ and thus complete the definition of the zigzag functions. 
Specifically: 

• G a method 

Exploring the possibility of g (k) to represent (f) {k ] and setting 

€>sl k) = GJQ^-\ (30) 

it is observed that Eqs. (21) and (30) are identical. Moreover, by invoking Eq. (7), the definition for G a 
given by Eq. (22) is obtained. 

• G Ia method 

When g\ k J , given by Eq. (29), is taken as a representation for (f {k) : , there results 

€>S^ = QfjG Ia - 1 (31) 

Substituting this expression into Eq. (7) yields yet another set of weighted-average constants 

G/-=lZ* W fiS < 32 ) 

h k = l 

and, consequently, with the use of Eqs. (13) and (23.2), another set of zigzag functions is obtained. 

Note that G a and G Ia represent two different weighted-average transverse shear properties of a 
laminate. In both cases, however, there exists dependence on the principal transverse shear moduli of each 
individual layer as well as their orientation angles 0^ k ^ . 

Alternatively, the homogeneous-limit strategy can be applied by examining the transverse-shear 
compliance coefficients instead of the stiffness coefficients. Hooke’s relations for the k-th material layer 
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in terms of the transformed compliance coefficients that relate the transverse shear strains to their 
conjugate stresses may be written as 



Therefore, the homogeneous-limit conditions for </)£. ^ and the compliance coefficients corresponding to 
the x a planes may be written 

if s£>^C a (34) 

Following the procedure outlined in Eqs. (28)-(32), two additional functions can be considered to 

i(k) • 

represent (p a ' , i.e., 

(«?> , , S%IC Ia - 1 ) (35) 

The two additional possible representations for the zigzag functions are 
• C a method 

Herein c { a k> is set to represent , i.e., 

= 1 (36) 

Substituting Eq. (36) into Eq. (7) yields 




(37) 


• C- Ia method 

Letting cfj to represent , there results 

€>^=s^ic Ia - 1 

Substituting Eq. (38) into Eq. (7) yields 


(38) 
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c z-=rf;* w 32 < 39 > 

where it is noted that C a and C Ia are two different weighted-average laminate constants for a given 
laminate. For the special case of a diagonal compliance/stiffness matrix representing all layers 
( Qu ^ = 0 ), such as in a cross-ply laminate whose principal material directions are aligned with the 
Cartesian plate coordinates, it can be readily verified that 

(C. , C, a ) = {G-l , G~') and ( C <*> . c«)=(g* . g?) (40) 

The four zigzag-function methods are summarized in Table 1. In Section 5, these methods are 
examined quantitatively by way of analytic solutions for simply supported sandwich plates. 


Table 1: Zigzag function definitions 


Homogeneous 
limit method 

Weighted-average 

constants 

Zigzag-function derivatives 

</> (k) 

Ta,z 

Interfacial 

displacements 

1 


\ 1 " /> <H ' 
^ n k= i y aa , 

-l 

g {k) =GjQ'y-\ 

(N 

II 

2 

G >« 

n k = 1 

g%=Q { "'G Ia -\ 

3 


( 1 " h (k) ' 

j'hs^y 

-l 

r— H 
1 

s§ 

to 

<0 

II 

V 

4 

c Ia 

r (*) - ?(*) / r _ l 

C Ia - 1 L la 1 


4.3 Modeling homogeneous plates with the full power of zigzag kinematics 

When modeling homogeneous plates or laminated plates with the same transverse shear stiffnesses on 
the planes of rectilinear orthotropy, the refined zigzag theory can be successfully applied by invoking the 
homogeneous limit property. To achieve accurate predictions for homogeneous plates without the use of 
shear correction factors, the differences in the layer transverse shear moduli can be made arbitrarily small, 
with such an approach taking the full advantage of zigzag kinematics. A judicious choice of small 
perturbations to describe nearly homogeneous transverse shear properties over the plate’s thickness can 
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lead to remarkably accurate predictions, achieving parabolic shear strain and stress distributions that are 
consistent with elasticity theory predictions for homogeneous plates. Moreover, the approach yields 
accurate nonlinear inplane strains and stresses when modeling thick plates. By contrast, the use of purely 
homogeneous transverse shear properties would result in the vanishing of the zigzag functions and would 
lead to the exclusively coarse theory, known as FSDT; where, for homogeneous plates, the kinematic 
assumptions give rise to constant through-the-thickness transverse shear strains and stresses and linear 
inplane strains and stresses, thus requiring a shear correction factor to compensate, in an average sense, 
for the lack of the parabolic transverse shear distributions. 

The approach is demonstrated by using the refined zigzag theory based on the G a method (see Table 
1). For simplicity and clarity of the discussion, a laminated plate whose principal axes of rectilinear 
orthotropy are coincident with theXj and x~, plate axes (thus =G[\\ Q 22 = G\t\ =0) is 
considered; further, cylindrical bending about the x 2 axis is considered so that the response quantities are 
functions of only the x l and z coordinates. Using Eq. (10.4), the refined zigzag theory yields the 
following expression for the transverse shear strain 



= h(x i) 


1 + ^)ZA> 

" M* i) 


(41) 


It is then assumed that the material layers have slightly different shear moduli from a constant value G n 
G®=G 13 (l + £<*>) (k = \,...,N) (42) 

where s (k) =SjU (k) 1 with /J (k) =0(1) denoting a dimensionless piecewise-constant function and s 
an arbitrarily small scalar. Substituting Eq. (42) into Eqs. (22) and using Eq.(21) yields 

G, = G l3 and $ k) - -SjLi {k) (43) 

Note that both y/ x (jc, ) and y x (x, ) have the same functional fonn in terms of x l ; also, for this nearly 

homogeneous case their magnitudes are respectively 0(s _1 ) and 0(1) . Therefore, the ratio of the two 
functions appearing in Eq. (41) is simply a large constant 


K(*i) 


OCT 1 ) 


(44) 


Accounting for Eqs. (43) and (44) yields the transverse shear strain of the form 

=ri(*i)(i-‘y' <t) ) (45) 

where c is a dimensionless constant and (1 — c// A) ) represents a piecewise constant thickness 
distribution of the transverse shear strain. 
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A parabolic transverse shear distribution is a good approximation for homogeneous plates that range 
from moderately thick to thin plates. Therefore, for the cylindrical bending problem considered herein a 
parabolic distribution of the transverse shear strain may be expressed as 


7\z — ri(*i) 



(46) 


where the r^Xj) function is determined from a given boundary value problem. 

The function // (A) in Eq. (45) can now be determined in such a way as to allow (1 -c// (A) ) to 
approximate a parabolic distribution given by Eq. (46) in an average sense. Thus, the average shear strain 
across the Mi layer resulting from Eq. (46) is set to correspond to (1 - c// <A> ) , i.e., 

1 f z w 

2 /><*> J -V.) 


1 -4L z=1 - 4- 1 ) +4) + Ww =1 _ ( 

h 2 3h 2 


from which 


(22 \ 
„(*)__?_ Z (*-l) Z W Z (k-\) Z (,k) 

M 3 c h 2 h 2 h 2 

V 

and, finally, 



(48) 


(49) 


where the constant 3c has been absorbed into an arbitrarily small scalar s 1 . 

A similar analysis of cylindrical bending about the A‘ 2 results in the expressions for the 
G 2 J transverse shear moduli 

G®=G 23 (i + £ W ) (50) 

where c (A) is the same as that given in Eq. (49). 

Equations (42), (49) and (50) provide the slightly perturbed transverse-shear-modulus distributions 
through the thickness, for the purpose of achieving a homogeneous plate in a limiting sense as the small 
parameter £ (k) diminishes to zero and the number of discretized layers N becomes sufficiently large. 
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4.4 Remarks on equilibrium equations, boundary conditions, and FEM approximations 

The plate equilibrium equations and consistent boundary conditions are derived in [30] from the 
virtual work principle, resulting in a system of seven second-order partial differential equilibrium 
equations in terms of the seven kinematic variables, thus constituting a 14 th -order theory. The number of 
kinematic variables within this theory is fixed at seven and does not depend on the number of discretized 
layers; hence the computational effort is practically the same for any number N . From the computational 
perspective of finite element approximations, -continuous kinematic functions can be used because the 
strain measures are represented by partial derivatives not exceeding first order. 


5. Example problems and results 


Analytic solutions for simply supported square laminates are used to examine the predictive capability 
of the four variants of the present Refined Zigzag Theory (RZT). The plate is defined on the domain 

Xj e x 2 e [0,a], z e [— /?, //] and is subjected to the sinusoidal transverse loading 

q = q 0 sin (xx x / tf)sin(/rx 2 la). The analytic solutions for uniaxial, cross-ply, and angle-ply laminates 

are derived in [30]. The first set of example problems examines the modeling of thin, moderately thick, 
and thick homogeneous orthotropic plates using the perturbed transverse-shear approach. Then, results for 
a highly heterogeneous and anisotropic angle-ply antisymmetric sandwich laminate are critically 
examined. The numerical and graphical results that follow include results of the four variants of RZT and 
several other theories that are used for comparison purposes; the results are labeled as: 


3D Elasticity: 


FSDT: 


Three-dimensional elasticity solutions using procedures developed by Pagano 
[36] for cross-ply multilayered plates and by Noor and Burton [37] for angle-ply 
antisymmetric laminates. 

First-order Shear Deformation Theory; shear correction factor k 2 = 5/6 . 


Reddy: 


Third-order equivalent single-layer theory by Reddy [9]. 


Zigzag (D): 


Di Sciuva theory [23]. 


RZT (G a ): 
RZT ( G a ; N): 

RZT (G la ): 
RZT ( C a ): 
RZT (C Ia ): 


Refined zigzag theory based on the G a method (refer to Table 1). 

Refined zigzag theory based on the G a method modeling a homogeneous plate 
with A discretized layers through the thickness (refer to Figures 3-5). 

Refined zigzag theory based on the G Ia method (refer to Table 1). 

Refined zigzag theory based on the C a method (refer to Table 1). 

Refined zigzag theory based on the C Ia method (refer to Table 1). 
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In Figures 3-8, the through-the-thickness distributions correspond to the normalized displacement and 
stresses, given as 


U\ =(10 4 D n /q Q a 4 ) u[ k) (0,a/2,z) , cr n = [(2h)' / q 0 a 2 ^ / 2,a /2,z) , and 

t 1z = (2h/q 0 a ) rjf (0, a/2,z). 


5.1 Homogeneous orthotropic plates 

The ability of RZT to model homogeneous orthotropic plates using the perturbed transverse-shear 
approach is demonstrated subsequently. The plate is made of a carbon-epoxy material with the 

mechanical properties given by Young’s moduli E[ k) = 1.579 x 10 2 GPa, E (k) = E (k) =9.584 GPa, Poisson 
ratios =0.32, vf i ) =QA9, and shear moduli =5.930 GPa, and G^ =3.227 GPa. For 

this homogeneous, orthotropic plate the principal material directions are aligned with the plate coordinate 
axes. 

The four variants of RZT, when applied to the carbon-epoxy (homogeneous and orthotropic) plates of 
various span-to-thickness ratios (a/2h=5, 10, and 100), showed close agreement, with RZT ( G a ) and 

RZT ( C Ia ) producing slightly superior results. For this reason, subsequent discussions concerning 
homogeneous plates focus exclusively on RZT ( G a ) that was originally introduced in [30]. The perturbed 

transverse-shear approach is assessed by fixing the small scalar at 5 = 10 5 while studying the plate’s 
response by increasing the number of the discretized layers, N . 

Results are shown for the maximum values of three response quantities in Table 2. These results 
include the central transverse displacement (deflection), which is averaged across the thickness, 

u a . ye [a / 2, a / 2) , the central inplane stress on the top surface, (al2,al 2,/?), and the edge 

transverse-shear stress at the midplane, rf. (0,tf /2,0). The results are normalized with respect to the 

corresponding solutions of three-dimensional elasticity theory. The response quantities are seen to 
converge from below as N is increased from 1 to 64. The special case of N =\ corresponds to FSDT 

( k~ =1), i.e., the zigzag functions are identically zero in this case. For a thin plate (0/2/7=100), the 
deflection and inplane stress match the exact solution very closely for all values of N , whereas the 
transverse shear stress converges to the exact solution of three-dimensional elasticity theory as N is 
increased. For moderately thick ( a/2h=\0 ) and thick (a/2h=5) plates, all quantities converge to slightly 
greater values that exceed those of three-dimensional elasticity. The converged results corresponding to 
RZT ( G a ; N=64) were also compared with Reddy’s third-order theory and both theories were found to be 

in close agreement. In both cases, the slightly over estimated results for thick plates may be attributed to 
the lack of transverse normal flexibility in these theories. 
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Table 2. Deflection, inplane stress, and transverse shear stress as a function of the number of 
discretization layers, N, for homogeneous orthotropic plates (a/2h= 5, 10, and 100); results correspond to 
RZT (G a , N) and are normalized with respect to three-dimensional elasticity solutions. 


N 

a/2h = 5 

a/2h= 1 0 

a/2h=\00 

uf 


4 



< 

u™ 


< 

1 

0.948 

0.870 

0.591 

0.976 

0.964 

0.565 

1.000 

1.000 

0.556 

4 

1.002 

0.972 

0.946 

0.999 

0.993 

0.925 

1.000 

1.000 

0.917 

16 

1.024 

1.015 

1.025 

1.009 

1.005 

1.003 

1.000 

1.000 

0.995 

64 

1.026 

1.017 

1.030 

1.009 

1.006 

1.008 

1.000 

1.000 

1.000 


Normalized through-the-thickness distributions for the u x , a n , and r lz response quantities, 

corresponding to a thick orthotropic plate (a/2h=5), are depicted in Figures 3-5. The special case of N=l 
shows a linear variation through the thickness that matches that of FSDT. As N is increased, a nonlinear 
distribution is achieved (Figure 3 (a)). A close comparison of the nonlinear, converged solution to the 
corresponding solutions of three-dimensional elasticity and Reddy’s third-order theory is shown in Figure 
3 (b). For this thick orthotropic plate, Reddy’s theory is perfectly applicable and its predictions are in 
very close agreement with RZT ( G a ; 7V=64). In Figures 4 and 5, similar comparisons are provided for the 

inplane stress, <t h , and transverse shear stress, Tj. , distributions. Of particular interest is the 

convergence of transverse shear stress; both RZT ( G a ; A=64) and Reddy’s parabolic stresses slightly 

over estimate the maximum value at the midplane; however, RZT ( G a \ N= 64) is slightly closer to the 

three-dimensional elasticity solution. An additional study of thinner homogeneous orthotropic plates 
reveals that all three theories achieve full agreement, demonstrating correct linear distributions for the 
inplane displacements and stresses and a parabolic variation of the transverse shear stress through the 
thickness (not shown). 

These results demonstrate that RZT is perfectly suited for predicting accurate parabolic transverse 
shear stresses, in a piecewise-constant fashion; hence no shear correction factors are required to model 
homogenous plates. Moreover, for thick plates, accurate nonlinear inplane displacements and stresses 
consistent with three-dimensional elasticity are obtained in a piecewise-linear fashion. Recall that, 
similarly, no shear correction factors are needed for heterogeneous composite and sandwich laminates, as 
discussed in [30], 
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5.2 Assessment of four zigzag-function methods for heterogeneous laminates 


This section focuses on the assessment of the four variants of RZT presented in Section 4. The 
selected laminate has a high degree of material heterogeneity and anisotropy and is the most challenging 
of the many laminates examined in [30]. The laminate is an angle-ply, antisymmetric, thick sandwich 
plate with the aspect ratio a/2h= 5. The five-layer sandwich has two stiff carbon-epoxy face sheets and a 
thick and very compliant PVC core. The nonnalized lamina-thickness distribution, h <k> /h, starting from the 
first layer, is given by (0.05/0.05/0.8/0.05/0.05); the corresponding lamina orientation angles are (30°/- 
45°/0°/45 o /-30 o ); the material distribution is given by (C/C/P/C/C), where the labels C and P correspond to 
the carbon-epoxy and PVC materials, respectively. The PVC material is modeled as an isotropic material 

with Young’s modulus £"=1.040 x 10" 1 GPa and Poisson ratio v=0.3. The mechanical properties of the 
carbon-epoxy plies (material C) are the same as those reported in Section 5.1. 

Comparisons of results for the center deflection, which is normalized with respect to the three- 
dimensional elasticity solution, are presented in Table 3. The results are obtained from FSDT, Reddy’s 
third-order theory, Di Sciuva’s zigzag theory, and the four variants of RZT. The results demonstrate that 
RZT ( G a ) and RZT (C la ) underestimate the deflection by less than 0.1%, compared to 12% by RZT 

( G Ia ) and RZT ( C a ). By comparison, FSDT underestimates the deflection by about 92%, and Reddy 

and Di Sciuva theories underestimate the deflection by about 62% and 10%, respectively. Note that Di 
Sciuva’s solution is only approximate, due to the presence of angle plies in this laminate (refer to [30] for 
further details). Also, FSDT’s deflection is expected to improve by employing lamination-appropriate 
shear correction factors. 


Table 3. Center deflection nonnalized with respect to three-dimensional elasticity solution for angle-ply, 
sandwich laminate (a/2h=5). 


Normalized center deflection, w Cellter(Theoiy) / w C enter(3D Elasticity) 

FSDT 

Reddy 

RZT (C„) 

RZT (G Ia ) 

Di Sciuva 

RZT ( G a ) 

Rzt ( C la ) 

0.075 

0.382 

0.879 

0.882 

0.902 

0.999 

0.999 
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Figure 3. Through- thickness distribution of inplane displacement for a homogeneous orthotropic plate 
(a/2h= 5). 




Figure 4. Through-thickness distribution of inplane stress for a homogeneous orthotropic plate ( a/2h=5 ). 
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Figure 5. Through-thickness distribution of transverse shear stress for a homogeneous orthotropic plate 
(a/2h= 5). 


Figures 6-8 demonstrate through-the-thickness distributions for the u x , cr n , and r lz quantities for the 
sandwich laminate where, for comparison, three-dimensional elasticity and Reddy theories are used. 
Figure 6 shows that RZT ( G a ) and RZT ( C Ia ) yield superior wj displacement predictions and, as 

expected, FSDT and Reddy theories are the least accurate. Figures 7 and 8 depict the inplane stress, <T n , 
and transverse shear stress, Tj. , distributions, respectively, where the comparisons are narrowed to the 
best performing zigzag theories, RZT ( G a ) and RZT ( C Ia ), and Reddy’s third-order theory. For both 
stress components, RZT ( G a ) and RZT (C Ia ) produce highly accurate results. By comparison, Reddy’s 

theory underestimates the maximum inplane stress by about 50% (Figure 7(b)), and over estimates the 
maximum transverse shear stress by a factor of 20 (Figure 8(b)); this latter result is typical for a higher- 
order theory applied to a sandwich analysis. Furthermore, a close examination of the transverse shear 
stress in the face sheets (Figure 8(a)) indicates that although both the RZT ( G a ) and RZT (C Ia ) 

predictions are nearly equally accurate, RZT ( C - Ia ) solutions in the face sheets are slightly superior. 

Integration of three-dimensional elasticity equilibrium equations, while invoking RZT’s inplane stresses, 
produces highly accurate, continuous through-the-thickness, transverse shear stresses (not shown; for 
details refer to [30]). 


24 




Thickness coordinate, z / h 



Figure 6. Through-thickness distribution of inplane displacement for an angle-ply, antisymmetric, 


sandwich laminate (a/2h= 5). 
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Figure 7. Through-thickness distribution of inplane stress for an angle-ply, antisymmetric, sandwich 
laminate (a/2h=5). 
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Figure 8. Through-thickness distribution of transverse-shear stress for an angle-ply, antisymmetric, 
sandwich laminate (a/2h=5). 


6. Conclusions 


The Refined Zigzag Theory (RZT) has been reformulated from a multi-scale perspective and produced 
superior results over a wide range of material systems and plate aspect ratios. Four sets of zigzag 
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functions, derived from the condition of limiting homogeneity of transverse-shear properties, provide 
viable modeling alternatives; however, only two sets of these functions guarantee superior predictions. 
For all material systems, there are no requirements for shear correction factors to yield accurate results. 
This variationally consistent theory, derived from the virtual work principle, requires simple C°- 
continuous kinematic approximations for developing computationally efficient finite elements. 

To model homogeneous plates effectively using the full power of zigzag kinematics, a multilayered 
modeling approach that employs infinitesimally perturbed transverse-shear stiffness properties has been 
demonstrated. The methodology pennits excellent predictions of all response quantities and does not 
increase the computational effort since the number of kinematic variables remains unchanged. 

Results of analytic solutions have been presented which reveal that RZT is a highly accurate theory 
over a wide range of span-to-thickness ratios and material systems, including very challenging sandwich 
plates that exhibit a high degree of transverse-shear flexibility, anisotropy, and heterogeneity. The theory 
is therefore ideally suited for large-scale finite element analyses and can be successfully used for 
designing high-performance homogeneous, laminated composite, and sandwich aerospace structures. 
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